#include "base.h"

int main(int argc,char **argv)
{


          //////////////////
          // Input   ///////
          //////////////////

          int Nr;
          char dir[1024];
	  int beta100;
	  char momentum_order[1024];
          int imass;


          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-imass",6)==0)  sscanf(argv[i],"-imass=%d",&imass);
            if(strncmp(argv[i],"-Nr",3)==0)  sscanf(argv[i],"-Nr=%d",&Nr);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100);
            if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);

	   
          }

		printf("====== imass = %d ===========\n", imass);
		printf("dir %s\n",dir);

   /////////////////////////
  //// Momenta set
  //////////////////////////////

  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  Momentum(momentum_order);
  // equivalent momenta
  Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
  Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);

 //////////////////////////////////
 // IO files
  ///////////////////////////////

  char file[1024];
	 

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

	//////////////////////////////////
	// Read Ds
       // I work with reals from now on
	/////////////////////////////

     printf(" Read Ds         \n");

     rO55 **DpPCr=malloc(sizeof(rO55*)*(Nmom));
     for (int imom=0; imom<Nmom;imom++) DpPCr[imom]=malloc(sizeof(rO55)*(Nclust+1));
	 

     complex *aux_clust=malloc(sizeof(complex)*(Nclust+1));
     char dir_r[1024];

    for (int imom=0; imom<Nmom; imom++)
     for (int irc=0; irc<2; irc++)
     {
          int *rvector=malloc(sizeof(int)*(4)); whichr_sort(irc, rvector);
          int r1,r2,r3,r4; r1=rvector[0]; r2=rvector[1]; r3=rvector[2]; r4=rvector[3];
      	  sprintf(dir_r,"%s/%d%d%d%d/",dir,r1,r2,r3,r4);
      	  for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
      	  {
                sprintf(file,"DpPC%d%d_m%d",iop1+1,iop2+1,imass);
		read_input_onemomenta_complex(dir_r,file,aux_clust,imom);
		for (int iclust=0; iclust<Nclust+1; iclust++) DpPCr[imom][iclust][irc][iop1][iop2]=aux_clust[iclust];
          }
    }



     printf(" r promedium             \n");

     O55 **DpPC=malloc(sizeof(O55*)*(Nmom));
     for (int imom=0; imom<Nmom;imom++) DpPC[imom]=malloc(sizeof(O55)*(Nclust+1));


     for (int imom=0; imom<Nmom; imom++) 
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
     {
	for (int iclust=0; iclust<Nclust+1; iclust++) 
		DpPC[imom][iclust][iop1][iop2]=(DpPCr[imom][iclust][0][iop1][iop2]+DpPCr[imom][iclust][1][iop1][iop2])/2.;
     }





        //////////////////////////////////////////////
        //     EQUIVALENT MOMENTA
        ///////////////////////////////////////////////
     printf(" PROMEDIUM OVER EQUIVALENT SINP2 \n");


     Equivalent_Momenta(SinP2);
	//for (int imom=0; imom<Nmom_eq; imom++) printf("%d %8.5f \n", imom, SinP2_eq[imom]);

     O55 **DpPC_eq=malloc(sizeof(O55*)*(Nmom_eq));
     for (int imom=0; imom<Nmom_eq;imom++) DpPC_eq[imom]=malloc(sizeof(O55)*(Nclust+1));


     double *aux_mom=malloc(sizeof(double)*Nmom);
     double *aux_eq_mom=malloc(sizeof(double)*Nmom_eq);
  
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
       for (int iclust=0; iclust<Nclust+1; iclust++)
       {
	 for (int imom=0; imom<Nmom; imom++) aux_mom[imom]= creal(DpPC[imom][iclust][iop1][iop2]);
	 Average_over_Equivalent_Momenta(SinP2,aux_mom,aux_eq_mom);
	 for (int imom=0; imom<Nmom_eq; imom++)  DpPC_eq[imom][iclust][iop1][iop2]=aux_eq_mom[imom];
       }
	 
     //////////////////////////////////////////
     /////// WRITE FILES
     //////////////////////////////////////////
     printf(" Write files        \n");

     char dir_qav[1024];
     sprintf(dir_qav,"%s/QAV/",dir);
     double *aux_clust_real=malloc(sizeof(double)*(Nclust+1));     
 
     for (int imom=0; imom<Nmom_eq; imom++)
     {
     // write .dat file with all D_ij and its error
     sprintf(file,"%s/p%03d/DpPC_m%d.dat",dir_qav,imom,imass);
     FILE *fout=fopen(file,"w");
     if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
     }
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
      {
                 for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]= creal(DpPC_eq[imom][iclust][iop1][iop2]);
                 fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],creal(DpPC_eq[imom][Nclust][iop1][iop2]),error_jack_real(aux_clust_real,Nclust));
      }
      fclose(fout);

      // write each D_ij in different file for all cluster
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
      {
                sprintf(file,"DpPC%d%d_m%d",iop1+1,iop2+1,imass);
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=creal(DpPC_eq[imom][iclust][iop1][iop2]);
                write_output_onemomenta_real(dir_qav,file,aux_clust_real,SinP2_eq,imom);
      }
     }


	// WRITE ALLP files (for plots)

      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
      {
	sprintf(file,"%s/D_ALLP/DpPC%d%d_m%d.dat",dir_qav,iop1+1,iop2+1,imass);

        FILE *fout=fopen(file,"w");
        if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
        }
        for (int imom=0; imom<Nmom_eq; imom++){for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]= creal(DpPC_eq[imom][iclust][iop1][iop2]);  fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],creal(DpPC_eq[imom][Nclust][iop1][iop2]),error_jack_real(aux_clust_real,Nclust));}
	fclose(fout);
      }
    return 0;
}
